##############################
##Programme Laden
#############################
 
### Verzeichnist mit Modulen
    MODULE.DIR <- "C:/Documents and Settings/Stefan_Hoderlein/Desktop/"

### Einbinden der Funktionen
    source( paste(MODULE.DIR, "Gdata1.R", sep="") )
     source( paste(MODULE.DIR, "gen.af.R", sep="") )
 
 
 #########################################################################
 ## Daten aufbereiten                                                     
 #########################################################################
 filename <- "C:/Documents and Settings/Stefan_Hoderlein/Desktop/newGroups40.csv"
 
 A1 <- read.csv(filename , sep=";", dec=",")

##########HH Char
Zett <- cbind(A1$agehd1,A1$washmach,A1$ncars,A1$ifshhtyp)

Zschl <- gen.af(af=Zett,no=1)

##########Control Function

Endodat <- cbind(A1$totexp,A1$ln.hhincome,Zschl)
Endopos <- Endodat[,-1]
H1 <- c(15,0.5,0.3)
#Vau <-rep(NA,length(A1[,1]))

Endomod <- lm(Endodat[,1]~ Endodat[,2] + Endodat[,3])

summary(Endomod)

Vau <-Endomod$res

#Vau<-rnorm(mean=0,sd=1,n=length(A1[,1]))

summary(Vau)

#summary(A1)
Data <- cbind(A1$year,A1$FOOD/A1$relFOOD*100,A1$totexp,A1$relFOOD/100,Zschl,Vau)

#Data <- cbind(A1$year,A1$FOOD/A1$relFOOD*100,A1$totexp,A1$relFOOD/100)

colnames(Data)<-c("Year","Food","TotExp","Fprice","Princomp","controlF")
summary(Data)

Preis <-Data[,4]

for (j in 1:length(Data[,1]))

    {
    
    for (i in 1:27)

        {
        if (Data[j,1]==(1973+i))    
        Data[j,4] <- mean(Preis[Data[,1]==(1973+i)])
        }
    
    }

    
#################Group data so that we have 9 different scenarios

Daten1 <- matrix(Data[(Data[,1]==1974)|(Data[,1]==1975)|(Data[,1]==1976)],ncol=6)
Daten2 <- matrix(Data[(Data[,1]==1977)|(Data[,1]==1978)|(Data[,1]==1979)],ncol=6)
Daten3 <- matrix(Data[(Data[,1]==1980)|(Data[,1]==1981)|(Data[,1]==1982)],ncol=6)
Daten4 <- matrix(Data[(Data[,1]==1983)|(Data[,1]==1984)|(Data[,1]==1985)],ncol=6)
Daten5 <- matrix(Data[(Data[,1]==1986)|(Data[,1]==1987)|(Data[,1]==1988)],ncol=6)
Daten6 <- matrix(Data[(Data[,1]==1989)|(Data[,1]==1990)|(Data[,1]==1991)],ncol=6)
Daten7 <- matrix(Data[(Data[,1]==1992)|(Data[,1]==1993)|(Data[,1]==1994)],ncol=6)
Daten8 <- matrix(Data[(Data[,1]==1995)|(Data[,1]==1996)|(Data[,1]==1997)],ncol=6)
Daten9 <- matrix(Data[(Data[,1]==1998)|(Data[,1]==1999)|(Data[,1]==2000)],ncol=6)


Daten1[,4] <- rep(mean(Daten1[,4][Daten1[,1]==1975]),length(Daten1[,1]))
Daten2[,4] <- rep(mean(Daten2[,4][Daten2[,1]==1978]),length(Daten2[,1]))
Daten3[,4] <- rep(mean(Daten3[,4][Daten3[,1]==1981]),length(Daten3[,1]))
Daten4[,4] <- rep(mean(Daten4[,4][Daten4[,1]==1984]),length(Daten4[,1]))
Daten5[,4] <- rep(mean(Daten5[,4][Daten5[,1]==1987]),length(Daten5[,1]))
Daten6[,4] <- rep(mean(Daten6[,4][Daten6[,1]==1990]),length(Daten6[,1]))
Daten7[,4] <- rep(mean(Daten7[,4][Daten7[,1]==1993]),length(Daten7[,1]))
Daten8[,4] <- rep(mean(Daten8[,4][Daten8[,1]==1996]),length(Daten8[,1]))
Daten9[,4] <- rep(mean(Daten9[,4][Daten9[,1]==1999]),length(Daten9[,1]))

summary(Daten1)

##################calculate price change

deltap1 <-  Daten2[1,4] - Daten1[1,4]
deltap2 <-  Daten3[1,4] - Daten2[1,4]                                                                                                                                                                                                                                                                                                                          
deltap3 <-  Daten4[1,4] - Daten3[1,4]
deltap4 <-  Daten5[1,4] - Daten4[1,4]
deltap5 <-  Daten6[1,4] - Daten5[1,4]
deltap6 <-  Daten7[1,4] - Daten6[1,4]                                                                                                                                                                                                                                                                                                                          
deltap7 <-  Daten8[1,4] - Daten7[1,4]
deltap8 <-  Daten9[1,4] - Daten8[1,4] 

print(c(deltap1,deltap2,deltap3,deltap4,deltap5,deltap6,deltap7,deltap8))

#########prepare the variables

Q2maldeltap2 <- Daten2[,2]*deltap2
Q3maldeltap2 <- Daten3[,2]*deltap2

###############This is the Left hand side

print(summary(Q2maldeltap2))

B1<- sort(Daten3[,3],index.return=T)
B2<-B1$ix

IndeX <- rep(NA,10)

for (s in 0:9)
    
    {IndeX[s+1] <- round(length(Daten3[,3])*(5+10*s)/100)}
    
Index <- B2[IndeX]
Quant3 <- Daten3[Index,3]
Quant3


################# End Quant t+1

B3<- sort(Daten2[,3],index.return=T)
B4<-B3$ix

IndeX <- rep(NA,10)

for (s in 0:9)
    
    {IndeX[s+1] <- round(length(Daten2[,3])*(5+10*s)/100)}
    
Index2 <- B4[IndeX]
Quant2 <- Daten2[Index2,3]
Quant2

############### End Quant t

Quant2 <- 0.5*(Quant3 + Quant2)



############These are points on Y2 vector


#Positions <- cbind(rep(Quant2,10),c(rep(Quant3[1],10),rep(Quant3[2],10),rep(Quant3[3],10),rep(Quant3[4],10),rep(Quant3[5],10),rep(Quant3[6],10),rep(Quant3[7],10),rep(Quant3[8],10),rep(Quant3[9],10),rep(Quant3[10],10)))

Quant4 <- Quant2 - 1.1



for (s in 1:9)

    {Quant4<-c(Quant4,Quant2 - 1.1 + 0.1*s)}
  
    Positions <- cbind(rep(Quant2,10),Quant4)
 

Diff <- Positions[,2] - Positions[,1]
Position <- cbind(Positions,Diff,rep(0,100),rep(-5,100))
Position
H1<-c(20,1.5,20)


#Y2minY1<-matrix(NA,ncol=10,nrow=10) 
#for (i in 1:10)
#{Y2 <- Quant3[i]
#    for (j in 1:10)
#       {Y2minY1[j,i] <- Y2 - Quant2[j]
#        Positions[,3] <- Y2 - Quant2[j]}
#    }

#print(summary(Y2minY1)
print(summary(Q2maldeltap2))

##############These are Differences between Y2 and Y1

#########Now the probs


P2minP1<-matrix(NA,nrow=100,ncol=4) 

#P2minP1<-matrix(NA,nrow=100,ncol=6)

for (k in 1:100)

   {    print(k)
        Ind1 <- sign(Position[k,3]-Q2maldeltap2)
        Ind1[Ind1<0] <- 0  
        Dat1 <- cbind(Ind1,Daten2[,c(3,5,6)])
        
        Erg1 <- Gdata1(data=Dat1,h=H1,x0=c(Position[k,1],0,0),noend=1) 
        Erg1 <- Erg1[1]
        
        ##########Bootstrap 
        #for (l in 1:101)
        
        #    {x<-c(1:length(Dat1[,1]))
        #    IndX<-sample(x,replace=T)
        #    Dat1boot<-Dat1[IndX,]
            
            
        Ind2 <- sign(Position[k,3]-Q3maldeltap2)
        Ind2[Ind2<0] <- 0 
        Dat2 <- cbind(Ind2,Daten3[,c(3,5,6)])
        
        Erg2 <- Gdata1(data=Dat2,h=H1,x0=c(Position[k,2],0,0),noend=1) 
        Erg2 <- Erg2[1]
       
        P2minP1[k,] <- c(Erg1,Erg2,max(Erg1-Erg2,0),min(Erg1,1-Erg2))   
   }
    
ERGEBNIS<-cbind(Position,P2minP1)
ERGEBNIS
matrix(ERGEBNIS[ERGEBNIS[,8]>0.01|ERGEBNIS[,9]>0.01],ncol=9)

plot(density(ERGEBNIS[,8][ERGEBNIS[,9]>0.01]))
